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Abstract 

This paper is a mathematical study of some aspects of the signaUing 
pathway leading to the activation of the transcription factor NFAT (nu- 
clear factor of activated T cells). Activation takes place by dephospho- 
rylation at multiple sites. This has been modelled by Salazar and Hofer 
using a large system of ordinary differential equations depending on many 
parameters. With the help of chemical reaction network theory we show 
that for any choice of the parameters this system has a unique station- 
ary solution for each value of the conserved quantity given by the total 
amount of NFAT and that all solutions converge to this stationary solu- 
tion at late times. The dephosphorylation is carried out by calcineurin, 
which in turn is activated by a rise in calcium concentration. We study 
the way in which the dynamics of the calcium concentration influences 
NFAT activation, an issue also considered by Salazar and Hofer with the 
help of a model arising from work of Somogyi and Stucki. Criteria are 
obtained for convergence to equilibrium of solutions of the model for the 
calcium concentration. 



1 Introduction 

The phenomena modelled mathematically in this paper are mechanisms which 
are part of the way the immune system works at the molecular level. For 
background on immunology the reader is referred to |T7] or [3T]. T cells are 
among the most important components of the immune system. They have the 
task of recognizing certain antigens and reacting appropriately. More precisely, a 
T cell recognizes a peptide (small protein) in combination with an MHC (major 
histocompatibility complex) molecule. The recognition takes place through a 
surface molecule, the T cell receptor. In what follows attention will be confined 
to T helper (Th) cells although some of the statements made may also apply 
to other types of T cells. In order for the T cell to be activated a second 
signal is also necessary. This comes from another surface molecule, CD28, which 
recognizes the molecules B7.1 and B7.2 on the antigen presenting cell carrying 
the peptide-MHC complex. The information about these recognition events 
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is propagated to the nucleus through various signahing pathways. The result 
is that the transcription factors NFAT, NFkB and AP-1 bind to the DNA, 
leading to the production of the cytokine IL-2 (interleukin 2). Much remains to 
be learned about these signalling pathways and mathematical modelling has a 
great potential to contribute to obtaining a better understanding of them. 

In the following attention will be concentrated on the part of the signalling 
network relating to NFAT (nuclear factor of activated T cells). It should be 
noted that although the abbreviation NFAT refers to T cells this transcription 
factor is important for signalling in many other types of cells. There are five 
different NFAT molecules and the one of relevance in what follows is that known 
as NFATc2 or NFATl. A model for NFAT signalling in T cells was introduced 
by Salazar and Hofer [22]. In fact they only deal with part of the pathway. An 
important stage in signalling is when there is a flow of calcium ions into the 
cytosol. During the activation of T cells this occurs when IP3 (inositol 1,4,5- 
trisphosphate) binds to receptors in the endoplasmic reticulum (ER), opening 
calcium channels and thus allowing calcium ions to flow down their concentra- 
tion gradient. This can be simulated experimentally by treating the cells with 
ionomycin, which leads to transport of calcium ions across membranes. The flrst 
step in the NFAT pathway included in the work of [35] and the first one to play 
a role in what follows, is this increase in the calcium concentration. The calcium 
binds to calcineurin, partly activating it. It also binds to calmodulin, which can 
then complete the activation of calcineurin. The activated calcineurin removes 
phosphate groups from NFAT, which is present in phosphorylated form in the 
cytosol of resting cells. The NFAT then undergoes a conformational change and 
moves to the nucleus where it can bind to DNA. The main model in [22] (which 
will be called the SH model in what follows) describes the dephosphorylation 
of NFAT and its transport between the cytosol and the nucleus. A subsidiary 
model describes the calcium influx. The aim of this paper is to obtain a deeper 
mathematical understanding of these models. 

The SH model is a system of 4N + 4 equations and contains + 4 param- 
eters, where N is the number of phosphorylation sites. The case of interest for 
NFAT is = 13 so that there are 56 equations and 134 parameters. Due to the 
large number of variables involved it might seem difficult to analyse the dynam- 
ical behaviour of general solutions of this system. Chemical reaction network 
theory (CRNT) [9] is a general tool for attacking this type of problem and it 
turns out to be very effective in this case. As we will show, one of its strongest 
theorems, the Deficiency Zero Theorem, can be applied to this system. The 
result is that for a given total amount of NFAT there is a unique stationary 
solution of the system and that every other solution converges to the stationary 
solution. 

The SH model describes the dephosphorylation of NFAT when the concen- 
tration of activated calcineurin is constant. In fact the calcium influx which 
leads to the activation of NFAT is a dynamical process. To assess the applica- 
bility of the SH model it is desirable to know whether the calcium concentration 
tends to a constant value at late times. This process is modelled in [22] by a 
two-dimensional dynamical system. It will be shown that for certain subsets 
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of the parameter space for this model the solutions do converge to a stationary 
solution. It is also shown that when this happens the long-time behaviour of the 
amounts of the different forms of NFAT occurring in the SH model is that they 
converge to the values they converge to in the SH model with an appropriate 
choice of parameters. 

The paper is organized as follows. Section [5] contains some basic material 
about chemical reaction network theory. The dynamical analysis of the SH 
model for NFAT phosphorylation with constant stimulation is in section [31 The 
dynamics of the calcium influx is investigated in section |31 The last section of 
the main text gives conclusions and an outlook. In an appendix the SH model 
is compared with a model for the NFAT signalling pathway defined in |10j . 

2 Chemical reaction network theory 

Chemical reaction network theory is a collection of methods for studying the 
dynamics of solutions of ordinary differential equations modelling systems of 
chemical reactions. Some concepts of this theory will now be reviewed. In 
CRNT the basic objects are finite sets S of species, C of complexes and TZ of 
reactions. The elements of C are formal linear combinations of elements of S 
with positive integer coefficients while the elements of TZ are ordered pairs of 
elements of C. The number of elements of S, C and TZ are denoted by to, n 
and r respectively. The set S consists of the substances taking part in the 
chemical reactions and in the example of the SH system it consists of 4(7V + 1) 
states of NFAT. The complexes are the combinations of species occurring on 
the left and right hand sides of the reactions. In the SH model each complex 
is just a single species. The reactions are ordered pairs of complexes, each 
representing the input and output of one reaction. In the case of the SH model 
they can be identified with ordered pairs of species. (Note that a reaction and 
the reverse reaction are counted separately.) Given a reaction network let Cs be 
the concentration of the species with index s. Consider a system of ordinary 
differential equations of the form 

^ = E riy^y'Ws-ys)- (i) 

{y.y')e-R. 

Here r{y,y') > are the reaction rates. They are assumed non-negative. The 
choice of these functions is often referred to as the kinetics. In this section only 
the most standard choice of kinetics will be considered. This is mass-action 
kinetics, where 

r{y,y') = kyy,cy. (2) 

Here kyy> are positive constants called the rate constants and — Yises "^s" ■ 
The equation ([1]) will be abbreviated to c = /(c). Here c is a vector of con- 
centrations Cs and so is a point of R"*. The positive and non- negative orthants 
are defined to be the sets of points of R™ whose coordinates are positive and 
non-negative, respectively. The quantity c is said to be positive (non-negative) 
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if it lies in the positive (non-negative) orthant. Because of its interpretation in 
terms of concentrations c should be non-negative in order to be of relevance for 
applications. 

The positive and non-negative orthants are invariant under the evolution 
defined by the ordinary differential equations of a chemical reaction network. 
The invariance of the non-negative orthant follows from that of the positive 
orthant by continuity. The invariance of the positive orthant is a consequence 
of a lemma which will now be proved. (Cf. Lemma II. 1 of '25! for a similar 
result.) 

Lemma 1 Consider a solution Cs{t) of ([T]) with the coefficients r{y,y') being 
given by mass-action kinetics. If Cs(to) > for some s and some time then 
Cs{t) > for all t > Iq for which the solution exists. 

Proof If the statement of the lemma is false then it can be assumed that Cs(ti) = 
for some ti > to- The time ti can be chosen so that Cs{t) > for all t < ti. 
The quantity Cg satisfies an equation of the form 

^ = -/_(c)c, + /+(c) (3) 

where /+ is non- negative. Since Cs is positive on the interval [to, ti) the inequal- 
ity 

|(logc.)>-/_(c) (4) 
holds. Integrating this equation and exponentiating gives 

Csih) > c,(to)exp (^-^ ' \U{cit))\d?j . (5) 

The integral in this expression is finite and so the inequality implies that Cg (ti) > 
0, a contradiction. This completes the proof of the lemma. 

The reaction network can be represented as a directed graph where the 
vertices are the complexes and edges represent reactions. The reaction network 
is said to be weakly reversible if whenever it is possible to link the complex y to 
the complex y' by a sequence of reactions it is also possible to link y' to y in the 
same way. In the case of the SH system the network is weakly reversible since 
in fact every reaction is reversible. The connected components of the reaction 
graph are called linkage classes and their number is denoted by Z. In the case of 
the SH model 1 = 1. An important object is the stoichiometric matrix TV. Its 
columns correspond to the reactions belonging to the network. The entries in 
a column are defined by the sums of coefficients of the different species in the 
complexes occurring in the reaction, with the coefficients on the left hand side 
being counted negatively and the coefficients on the right hand side positively. 
In other words, these are the net number of molecules of each species produced 
in the reaction. It is an m x r matrix. The cosets of the form c -f im TV are 
called stoichiometric compatibility classes and are invariant under the flow of 
the system. The intersection of a stoichiometric compatibility class with the 
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non-negative orthant is called a reaction simplex and is also invariant under 
the flow by Lemma 1. The rank of N (i.e. the dimension of the stoichiometric 
compatibility classes) is denoted by s. The deficiency of the network is defined 
hy S = n — I — s. The following is part of the Deficiency Zero Theorem which 
was first proved in [13], [TS] and 0- 

Theorem 1 Let c = /(c) be the system of ordinary differential equations de- 
fined by a chemical reaction network by means of mass-action kinetics. If the 
network is weakly reversible and of deficiency zero then there is a unique posi- 
tive stationary solution in each stoichiometric compatibility class. The solution 
is asymptotically stable within its class. 

An important part of the proof of Theorem 1 is to show that there exists 
a Lyapunov function L{c) which is non-increasing along solutions and strictly 
decreasing along all positive solutions except for the stationary solution. This 
means that the stationary solution is the only possible positive w-limit point of 
a positive solution. The function / can be written in the form Yg where Y is 
called the complex matrix. Its columns are in one to one correspondence with 
the complexes and the entries in the column corresponding to the complex y are 
the components j/s. In the case of the SH system Y is just the identity. If c* 
is a stationary solution then /(c,) = 0. If in addition .g(c,) = the stationary 
solution is called complex balanced. Evidently any stationary solution of the SH 
system is complex balanced. In fact it is a consequence of the proof of Theorem 
1 that for a system satisfying the assumptions of that theorem the intersection 
of the kernel of Y with the image of N is {0}, so that any stationary solution 
of a system of that type is complex balanced. 

3 The model for NFAT phosphorylation 

The basic variables in the model of 22\ are amounts of different forms of NFAT. 
Using amounts rather than concentrations avoids introducing extra factors of 
the ratio of the volumes of the two compartments. This is just a matter of 
mathematical convenience. Phosphate groups can be attached to this molecule 
at up to = 13 sites. The index n will be used to denote the number of 
phosphate groups and runs from zero to N. It is assumed that the phosphate 
groups are bound to the sites in a certain order and are removed in the reverse 
order. There are thus A^ phosphorylation states in total. Each of these has 
an active and an inactive form. Each of them occurs in the cytosol and in 
the nucleus. This gives a total of 4(A^ + 1) variables decribing the amounts of 
the different substances. The processes of attaching a phosphate group and the 
conformational change between the active and inactive forms are reversible. It is 
necessary to prescribe 6 A^ -I- 2 rate constants to describe the reactions in a given 
compartment. However the rate constants describing the transitions between 
active and inactive forms are chosen to be the same in both compartments. Rate 
constants are also required to describe the transport processes between the two 
compartments - the active form is transported into the nucleus and the inactive 
form out of the nucleus. It is assumed that there is just one rate constant 
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for each of these two processes, independent of the phosphorylation state. A 
diagram of this reaction network can be found in |22j, Fig. 1. Assuming mass- 
action kinetics leads to a system of ordinary differential equations for the time 
evolution of the amounts of the different substances. 

The unknowns are as follows. The amount of active NFAT in the cytoplasm 
with n phosphorylated residues is denoted by a„, n = 0, 1, . . . , iV. The amount 
of the corresponding inactive form is denoted by i„. The amounts of these 
substances in the nucleus are denoted by An and /„. All these quantities are 
supposed non-negative. Mass-action kinetics is assumed. For 1 < n < — 1 
the dynamical equations for amounts in the nucleus are 

dAn 

— r— = Kn-lAn-l — KnAn + C„A„_|_i — C„_iA„ 

at 

+l-In-l+An + dan, (6) 

+l + An - I- In - fin (7) 

with rate constants Ar„, K'^, C„, C^, /+, /~, d and /. Evolution equations for the 
cases n = and n = N can be obtained by taking equations formally identical 
to those above with the conventions that K^i = C_i = K'_i = C'_i = and 
that Cn = C'j^ = Kn = K'^ — 0. Analogous evolution equations for the 
quantities a„ and i„ can be obtained as follows. Replace An and by a„ and 
in respectively everywhere except in the last term of each equation. Reverse the 
sign in the last term of each equation. Replace Cn, C'n, Kn, K'n by c™, c'n, kn 
and k'n- It is stated in [22] that the transport processes are much slower than the 
reactions within each compartment. On a heuristic level this can be imported 
into the mathematics by assuming that the coefficients d and / are very small. 
It may be hoped that solutions of the full system can be approximated by 
solutions of the system obtained by setting d = / = 0. In the latter system the 
equations describing amounts in the nucleus and the cytoplasm decouple. For 
this reason it is referred to in what follows as the decoupled system. To analyse 
this system it is enough to analyse the subsystems describing the dynamics in 
each compartment. 

The SH system satisfies n = 4:{N + 1) and 1^1. In order to show that 
Theorem 1 applies it therefore suffices to show that the rank of N is 4N + 3. 
This rank is equal to 4{N + 1) minus the dimension of the kernel of (N)'^ . 
The conditions for a vector to lie in this kernel are easy to analyse. Some of 
the conditions imply that all the components of the vector corresponding to 
amounts of substances in the cytosol are equal and that those corresponding to 
amounts of substances in the nucleus are equal. Finally elements corresponding 
to the amounts of the same substance in both compartments are equal. Thus the 
kernel is one-dimensional and the rank of the stoichiometric matrix is AN + 3. 
It follows that the deficiency is zero and Theorem 1 applies to the SH system. 

Theorem 1 leaves open the question whether under its hypotheses every 
solution converges to a stationary solution. It will be shown below that this is 
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the case for the SH system. Note first that the total amount of NFAT, which 
is the sum of all the variables in the SH system, is a conserved quantity. Thus 
each reaction simplex is compact. It follows that all solutions exist globally in 
time towards the future and the w-limit set of any positive solution is connected. 
Combining this with the remark about w-limit points made previously shows 
that if a solution does not converge to the stationary solution its w-limit set 
must be contained in the boundary of the positive orthant. 

Suppose now that a positive solution c has an w-limit point in the bound- 
ary of the positive orthant. There is a solution of the system, say Coo, passing 
through that w-limit point. The range of Coo is contained entirely in the bound- 
ary. By Lemma 1 the number of non-zero components of Ccc can never decrease. 
It might a priori increase but in any case it will be constant after a finite time. 
Thus when considering late-time behaviour it may be assumed without loss of 
generality than (coo)s(i) is non-zero for 1 < s < fc and identically zero for 
k + 1 < s < ra. Here < fc < to. In fact fc > since the sum of the vari- 
ables (coo)s, which is the total amount of NFAT in the cell, is conserved. If 
fc -I- 1 < s < TO then (coo)s = 0. There are no negative contributions to (coo)s 
since reactions having species s on their left hand side are not active. If there 
is a link from species s to species s' in the reaction network with non-zero con- 
centration then a positive contribution to (coo)s results. It follows that if s' is 
adjacent to s in the network then fc + 1 < s' < to. Since the reaction graph 
is connected this implies that all Cs vanish identically, a contradiction. Thus it 
has been proved that there can be no w-limit points on the boundary and the 
following result is obtained: 

Theorem 2 Let c = /(c) be the system of Salazar and Hofer. There is a 
unique stationary solution c* in each stoichiometric compatibility class and each 
positive solution converges to a stationary solution as t — oo. 

It has been conjectured that under the hypotheses of Theorem 1 every solu- 
tion converges to a stationary solution as t — ?■ oo. This is known as the global 
attractor conjecture 0] and has recently been proved in the case that there is 
only one linkage class by Anderson [2 . Theorem 2 could be deduced from the 
result of [2] but it has been shown here that there is a much easier proof in this 
relatively simple case. 

There are many different ways in which multiple phosphorylation can be 
organized and this can give rise to many systems related to the SH model. The 
phosphorylation is said to be processive if an enzyme which binds its substrate 
once phosphorylates several sites before dissociating. It is said to be distributive 
if only one site per binding event is phosphorylated. One type of distributive 
phosphorylation is sequential phosphorylation, where the sites are phosphory- 
lated in a particular order and dephosphorylated in the reverse order - this is the 
case in the SH model. There is also a cyclic variant where dephosphorylation 
takes place in the same order as phosphorylation. It is also possible to consider 
phosphorylation in a random order or mixtures of the mechanisms introduced. 
An extensive discussion of the possibilities and examples of biological systems 
where they occur can be found in [23l . For many of these systems an analogue 



7 



of Theorem 2 holds and can be proved in a similar way. For the only properties 
required for the proof are as follows: 

• each complex consists of one species 

• there is only one linkage class 

• the network is weakly reversible 

Thus the analogues of the SH model with cyclic or random phosphorylation both 
have the property that there is a unique stationary solution in each stoichio- 
metric compatibility class and that any other solution converges to a stationary 
solution. 

There is another dynamical system related to the modelling of T cell acti- 
vation which has deficiency zero and can thus be shown to have the property 
that any solution converges to a stationary solution and that there is only one 
stationary solution in any stoichiometric compatibility class. This is the kinetic 
proofreading model of McKeithan [TB] for antigen recognition by the T cell re- 
ceptor and its dynamics was analysed mathematically by Sontag [25]. In that 
paper the analogue of Theorem 2 is proved for McKeithan's model. The key 
observation is that the deficiency of the network is zero so that Theorem 1 ap- 
plies. From there it is possible to obtain the analogue of Theorem 2 for that 
system in a way very similar to what has just been done for the SH model. 

In the SH model each phosphorylation or dephosphorylation is modelled as a 
single reaction and the details of the interaction with the enzyme which catalyses 
the process are not included. Suppose that instead the enzyme is incorporated 
in the standard Michaelis-Menten way [18]. This means that the reactions 
describing the formation of a complex of the substrate with the enzyme and 
the dissociation of the complex to give either enzyme and substrate or enzyme 
and product are included, using mass-action kinetics. This is what is referred 
to as Michaelis-Menten via mass action (MMvma) in [TT] and is different from 
using an effective Michaelis-Menten kinetics for a single reaction. The analogue 
of one of the parts of the decoupled SH model with the simple mass-action 
kinetics replaced by MMvma kinetics is similar to what is called a multiple 
futile cycle in [26]. In that case there is only one kinase which catalyzes all 
phosphorylations and one phosphatase which catalyzes all dephosphorylations. 
This is slightly different from the situation in [11], where there is a different 
enzyme for each reaction. In |26| upper and lower bounds for the number of 
stationary solutions of a system of this type are obtained. It follows from these 
that while in the case = 1 there is only one stationary solution there are at 
least three stationary solutions for = 2 and at least thirteen for = 13 for 
suitable choices of the parameters of the system. This corresponds to the case 
where the total concentrations of the enzymes are small compared to the total 
concentration of the substrates. The number of stationary solutions is never 
greater than 2A^ — 1, whatever the parameters. If the total concentrations of 
the enzymes are sufficiently large compared to the total concentrations of the 
substrates then there is at most one stationary solution. 
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The equations arising in the decoupled system can be analysed in the same 
way as the full system and the analogue of Theorem 2 holds in that case. For 
the decoupled system more can be done and the stationary solutions can be 
calculated explicitly, as was shown in (35] . They are obtained by setting the sum 
of certain pairs of terms to zero. In the terminology of CRNT these stationary 
solutions are detailed balanced. Whether this gives the most general stationary 
solutions of the decoupled system is not discussed in [22] but it follows from 
the analogue of Theorem 2 that they are. To be concrete the system describing 
concentrations in the cytosol will be considered. To obtain the class of solutions 
found in [35] it is assumed that the first two terms on the second line of the 
evolution equation for a„ cancel. This also gives a similar cancellation in the 

evolution equation for i„. The condition for this is that ^ = L„ where L„ — 
Next it is assumed that the second and third terms in the evolution equations 
for a„ cancel for < n < iV — 1, giving 

^71+1 

Then the first and fourth terms cancel except in the case n = 0. In fact, if 
the equations involving L„ are satisfied and the equations (|S]) are satisfied for 
all n < — 1 then the conditions for a stationary solution of the part of the 
decoupled system describing concentrations in the cytosol is satisfied. Note that 
these detailed balance conditions cannot be satisfied by a stationary solution of 
the full system with non-zero coefficients d and /. They may, however, be 
approximately satisfied when d and / are small. For these stationary solutions 
of the decoupled system the fraction of the NFAT in the cytosol which is in the 
active state can be computed. It is given by 



AT 1 + fn""^^ 



(9) 



In order to have a better understanding of the system it is useful to consider 
a special case with a reduced number of parameters. This is obtained in the 
following way. The coefficients /c„, fc^, c„ and are taken to be independent of 
n and denoted by k, fc', c and c' respectively. It is also assumed that L„ = LqA". 
In this case the expression for (/> becomes 



v^JV /fey 



E;=o(i+ioA")(l)" 



L. 



.fef+l_i 



(10) 



One of the results of [H] is that for large N the function (j) resembles a Hill 
function 4'h{c) = a+c'^ fo'" ^ constant A and an exponent N. In what sense 
does this resemblance hold? If N is allowed to tend to infinity for fixed A then 
the Hill function tends pointwise almost everywhere to a translated Heaviside 
function which is zero for c < 1 and one for c > 1 . From the point of view of the 
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applications this gives rise to a switch behaviour for large N. For values of the 
control parameter c smaller than a threshold almost all the NFAT in the cytosol 
is in the inactive form while for values larger than the threshold almost all the 
NFAT is in the active form. The limiting behaviour of the function PU)) as N 
tends to infinity depends on the assumptions made about the other parameters 
present. If the other parameters are fixed then what is obtained in the limit 
does contain a threshold but is not a switch. The amount of activated NFAT 
is very small below the threshold but increases gradually above the threshold. 
This should be compared with the discussion in [T2] . Consider first the effect of 
varying c while keeping the other parameters fixed, (p (1 + Lo)~^ asc— >'00 
and — >• (1 + LqA^)"^ as c — >■ 0. The quantity A is assumed to be greater 
than one in [22^ and so if N is large the value of (f> at zero is close to zero. 
Consider next what happens if N tends to infinity for fixed values of the other 
parameters. In the region where c > Xk the inequality ^ < 1 holds. Then the 
limit of as oo is given by 



(11) 



c — Afc + Lo(c — k) 

When Afc tends to c the function (ioo tends to zero. Let c = c — Afc. Then 



(l + Lo)5 + Lofc(A- 1) 



(12) 



for c > Afc. In the region where c < Afc we get (f)^ = 0. Thus 0oo is a truncated 
translated Hill function with exponent one as mentioned in [T2]. There is also 
another interesting way of passing to the limit N —i' oo which is more closely 
related to what is done in [22]. To see this it is convenient to introduce the 

_2_ N_ 

variable fj, — Lq X. Inverting this gives Lq = (j) ^ . Substituting this into the 
expression for </> gives 



(13) 



Denote the limit of this function as 0oo- Assume that j is a fixed number 



which is less than one. When c > \/ — k the function is equal to one while 

when c < y^fc it is equal to zero. Thus i^oo is a translated Heaviside function 
and this limit behaves in a similar way to the limit of a Hill function when the 
exponent is allowed to tend to infinity while all other parameters are fixed. It 
is interesting to note that the threshold in occurs in a different place from 
the threshold in <j)oo- In the example plotted in Fig. 2(b) of '22' the parameters 
are chosen as A = 10 and /i = 1. Then the threshold for 0oo is at about 
3.3. At the threshold value of c the function </> has exactly the value one half. 
One key property implemented by this choice of parameters is that the reaction 
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constants are such that there is a strong tendency for weakly phosphorylated 
NFAT to change from the inactive to the active conformation and for highly 
phosphorylated NFAT to change from the active to the inactive one. 

The vaUdity of the above explicit formulae is restricted to the decoupled 
system but they do provide some information about the full system where d 
and / are non-zero. When all other parameters are fixed there is a unique 
stationary solution in a given reaction simplex for each choice of non-negative 
values of d and /. It can easily be shown, using the compactness of the reaction 
simplex, that the stationary solution depends continuously on the parameters. 
Since its position is known explicitly when d = / = its position is known 
approximately when these two parameters are small. 

For any solution of the SH system let (f> be the proportion of NFAT in the 
cytosol which is in the active state and ip the fraction of NFAT in the nucleus 
which is in the inactive state. Let Z be the fraction of NFAT which is in the 
nucleus. Let Vi and V2 be the volume of cytosol and nucleus. Then by the 
conservation of the total amount of NFAT the relation 



fZi;V2 = d(l - Z)cl)Vi (14) 
holds for any stationary solution and so 

^-d^v, + m,- ^^^> 

If / and d small then we can obtain approximate expressions for (j) f^nd ip as 
functions of the reaction constants. Thus a corresponding approximation is 
obtained for Z. 



4 Modelling the calcium influx 

The model presented in the last section gives a description of how the propor- 
tion Z of active NFAT in the nucleus depends on the parameters describing 
the state of the cell. An idealized experiment would then consist in modifying 
the rate constants c and C by stimulating the cell and measuring the resulting 
change in Z. In real experiments things are more complicated. Consider now 
the experiments carried out in ,19 . There, among other things, the following 
type of experiment was done. A population of T cells was treated with differ- 
ent concentrations of ionomycin and the production of IL-2 by these cells was 
measured. To get production of IL-2 the second signal is also required and it is 
produced artificially by stimulating the cells with PMA (phorbol 12-myristate 
13-acetate). The result for a given level of stimulation is a curve describing 
the number of cells producing different amounts of IL-2. If all the cells were 
identical this curve would reduce to a Dirac S but in reality this is smeared out 
to a smooth curve by the natural variability of the cells. When the stimulus is 
varied the following phenomenon is observed. The curve in general consists of 
two peaks separated by a region of production rates where there are very few 
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cells. As the stimulus is varied the height of one of the peaks grows while that 
of the other shrinks. The value at which the maximum of a peak is attained 
is essentially unchanged. The interpretation of this result is that a given cell 
either produces almost no IL-2 or produces IL-2 at a close to maximal rate. It 
exhibits a switch behaviour. This corresponds to the switch-like dependence of 
(j) on c found in the model. 

In the analysis of these results it is supposed that the number of cells where 
IL-2 production is switched on closely reflects the number of cells where Z is 
almost one. On the other hand the idea that different levels of stimulation can be 
modelled by taking different constant values of c and C may be oversimplified 
and this issue will now be looked at in more detail. The idealized situation 
would be to set a constant level of calcineurin activity by setting a constant 
calcium concentration. In fact it is observed experimentally that when a calcium 
influx is caused by stimulating the TCR and CD28, or by treating the cell with 
ionomycin and PMA, the concentration of calcium often displays oscillations. 
It is also possible that these oscillations, rather than being an unimportant 
side effect, encode information. A rise in calcium concentration does not only 
affect the NFAT signalling pathway but also other characteristics of the cell. 
This raises the question of how one signal can control several outputs. One 
possibility is that the information is encoded in the time dependence of the 
calcium concentrations and that the oscillations have an important role to play 
in this. Some experimental results on calcium oscillations in liver cells and 
T cells are described in [24] and [5] respectively. For an extensive review of 
modelling of calcium dynamics and signalling see [7] . 

The constants c and C represent the concentration of active calcineurin 
and so in modelling the effects of calcium influx they should be replaced by 
functions of time. How does the concentration of active calcineurin depend on 
the concentration of calcium? In ^22] this is modelled by an equation of the 
form 

dz 

— = a{zo - z)y - bz (16) 
at 

where z is the concentration of active calcineurin, y is the concentration of 
calcium in the cytosol and a and b are constants. It remains to model the 
dependence of the calcium concentration on time, possibly including oscillations. 
In |22j this is done by means of a system of two ordinary differential equations. 
These equations contain the concentration of IP3. With a view to modelling 
stimulation by ionomycin it will be assumed that the concentration of IP3 is 
constant. Then the model for the calcium concentration is schematically of the 
following form: 

X = p[-a{x -y)+ /3y- \f{y){x - y)], (17) 
y = a{x - y) - (3y + Xf{y){x -y)+'-f-Sy. (18) 

Here x is the concentration of calcium in the lumen of the endoplasmic reticu- 
lum. The quantities a, /?, 7, 5, A and p are positive constants and / is a positive 
function describing the response of the IP3 receptor to the calcium concentra- 
tion. In [22] it is chosen to be a Hill function with exponent two. This system is 
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essentially a special case of one introduced by Somogyi and Stucki [53] . The only 
reason it is not included is that the model in [24] has p = 1. Thus it implicitly 
assumes that the cytosol and the lumen of the endoplasmic reticulum have the 
same volume. In [23] another variant was also considered where the fact that y 
is much smaller than x is used to replace x — y hy x wherever it occurs in the 
above system. The above system will be called the modified Somogyi-Stucki 
model while the other variant will be called the unmodified Somogyi-Stucki 
model. Note that in the unmodified system a multiplicative constant like p can 
be removed by rescaling the variable x and some of the parameters. It is re- 
marked in [24] that taking f(y) — and a = in the unmodified model gives 
a system which is equivalent to the well-known Brusselator [20] . 

Now some remarks will be made on the dynamics of solutions of the system 
(fT7)) - p8)) . Note first that the only terms on the right hand side of the evolu- 
tion equation for a given unknown which are negative contain a factor of that 
unknown. Thus it can be shown as in the proof of Lemma 1 that a solution 
which starts positive remains positive. Taking a linear combination of the two 
evolution equations shows that any stationary solution (a;*,y*) satisfies y* = j. 
Substituting this back in the equation x = gives 

7[a + /? + A/(yO] 

S[a + Xf{y.)] ■ ^''^ 

Thus, in particular the system has exactly one stationary solution for any choice 
of the parameters and the function /. Linearizing the system about the station- 
ary solution gives 

— = p{[~a - Xf{y.)]S: + [a + /3 - Xf'iy.){x, - y,) + A/(j/,)]y}, (20) 

^ = [a + A/(y,)]i + [-a - /3 + Xf'{y.)ix, - y,) - Xf{y,) - S]y. (21) 

The determinant of the linearization at (a;*,j/*) is 6p{a + A/(j/*)) > 0. Thus 
the linearized stability of the stationary point is determined by the sign of the 
trace of the linearization which is given by 

-il + p)a~ 13-6^(1 + p)Xf{y,) + Xf'{y,)ix, - y,). (22) 

The stationary solution is a hyperbolic source if and only if 

p<{a + Xf{y,))-'[-a - fi - 5 - Xf{y,) + A/'(y,)(x, - y,)]. (23) 

Since and do not depend on p this shows that the region of instability is 
non-empty and is bounded by a smooth hypersurface in parameter space. 

2 

In the model of [22 the nonlinearity is given by f{y) = ^ . It follows that 
fiy) - and so 

~{i + p)f{y*) + .ny*)i^*-y*) 
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^ '^^U + vl^ {A + vl) 5[a{A + yl) + \yl)] 

^ {-(g + A)(l + p)yl + A[Al + P)a + W]}yl , . 

{A + yl)[a{A + yl) + Xyl] ' ^ > 

This means in particular that in this case a necessary condition for the positivity 
of the trace is 

/7X2 A[2^-£+p)a] 

\~5) < (a + A)(l+p) ■ ^^^^ 

Now it wiU be proved that every solution of the Somogyi-Stucki model is 
bounded. Note first that 

^ (x + py) = p(7 - 5y) . (26) 

Let Lc be the part of the line x + py — C for which y > j/5. No solution can 
cross it in the direction of increasing y. Let L' be the line y = ^jqrg^rjqrj- It 
cannot be crossed in the direction of decreasing y. Let H be the zero set of the 
right hand side of ([T7)) . Its equation can be written as 



y 



a + Xf{y) 



(27) 



This curve starts at the origin and as x tends to infinity along it y tends to 
infinity. Let 



-4 



1 + - 



A/(?) 



(28) 



Then the curve crosses the line y = j/6 when x = M. Define a region Rc to be 
that bounded by the y-axis, the line Lq, the vertical line joining the endpoint 
of Lc with L' and the part of L' between its intersection with that vertical 
line and the y-axis. For C > M + ^ this region is invariant under the flow. 
Hence any solution which starts in Rq remains there and is bounded. Since 
the union of the allowed Rc cover the region y > it follows that any 

solution which tended to infinity would have to lie in the strip y < ^jqrg^rjipx ■ 
Hence for a solution which is unbounded y > M implies y < 0, a contradiction, 
and all solutions are bounded. It follows from Poincare-Bendixson theory that 
unless the trace of the linearization at the stationary solution vanishes each 
solution either converges to a stationary solution or to a periodic solution. The 
condition on the trace is required so as to rule out possible orbits homoclinic to 
the stationary solution. 

For certain ranges of the parameters it can be shown that there are no 
periodic solutions. To see this note first that any periodic solution must lie in 
the region y > 2J+^a+3+A- '^^ '^^^ neither lie entirely in the complement of 
that region, due to the monotonicity of y there, nor leave that region. It must lie 
in some Rc and in fact it must do so for some C < M + ^. For the boundary of 
Rc with any larger value of C can only be crossed in one direction. In particular 
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for this solution x can never exceed M. Consider the divergence of the vector 
field defining the dynamical system. If the divergence is negative on a region 
containing the solution being considered then a contradiction is obtained. This 
divergence is given by replacing and j/* in (|22p by x and y respectively and 
can be written in the form 

-[{p+l)a + /3 + 5+{p+ l)Xf{y) + Xf'{y)y] + Xf'{y)x. (29) 

The last term is the only positive one and it can be bounded above by 3%/3AM 
when X < C. Hence if 

5^^^<(p+l)a + /3 + <5 (30) 

the divergence is negative on Re- It follows that a sufficient condition for the 
absence of periodic solutions is that 



3V3A7 / P 



8A-2S V a + A/(i) 



<{p + l)a + l3 + S. (31) 



When there are no periodic solutions it can be concluded that every solution 
converges to the stationary solution. 

In [22] numerical simulations were done for the above model with certain 
values of the parameters. These are: 

a = 0.001, 7 = 0.38, S^l, A = 0.16, p = 10, A = 0.25. (32) 

In two different simulations /3 was chosen to be 0.1 and 1. In both these cases the 
quantity ()22p is negative and so the stationary solution is a hyperbolic sink. For 
(3 — OA the inequality is violated and so this necessary condition suffices 
to verify that the stationary solution is a sink in that case. For /3 sufficiently 
large with the other parameters as above the trace is positive. The value of 
[3 for which the trace vanishes is about 8.2. For = 0.1 the inequality (|31|) 
is satisfied and the criterion for the absence of periodic solutions applies. For 
(3 ~ 1 it docs not. 

Suppose that for certain values of the parameters in (jl7p - (|18l) and certain ini- 
tial data the solution converges exponentially to the stationary solution (a;,, y^) 
as t — > 00. In other words, assume there are positive constants rj and R such 
that 

\x{t)-x,\ + \y{t)-y,\<Re-''*. (33) 



The equation (|16p can be rewritten as 

dz 

— ^ -(b + ay)z + azoy. (34) 

Define 

azay^ 

z* = 7- ■ (35) 

+ ay. 
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Then it follows from the variation of constants formula that 

zit) = ^(io)e-''('-*«'e-"^'>(^)''^ + re-''(*-^)e-"/X")'*"azo2/(s)ds (36) 

J to 

and hence that z{t) = + . . . where the remainder decays exponentially as 
t — >■ cx). 

What can be said about the asymptotic behaviour of solutions of the system 
obtained when the SH model is generalized by replacing the constants c and 
C by functions of t which converge to constants c* and C* respectively as t — ^ 
oo? Denote by X{t) a solution of the system of ODE obtained from the SH 
model by allowing the coefficients to be time dependent. Denote this system 
by X = F{X). Suppose that these coefficients converge to positive limits as 
t —i' oo and that their time derivatives tend to zero. If t„ is a sequence tending 
to infinity as n — c» consider the translates Xn{t) = X(t + tn). Because of 
the compactness of the reaction simplex these are uniformly bounded. Using 
this in the evolution equations shows that the time derivatives X'^{t) are also 
uniformly bounded. Differentiating the equation shows that the sequence X'^{t) 
is uniformly bounded. By the Arzela-Ascoli theorem there is a subsequence 
which converges together with its first derivatives uniformly on compact subsets 
to a limit Xoo{t). The functions Xn satisfy 

X^{t) = F{X^{t)) - F{X{t + tn)) = F^{Xn{t)) + .... (37) 

Here F^o is obtained by replacing the coefficients in F by their limiting values 
and the remainder term converges to zero as n — > oo. Hence 

XUt) = Foo{Xoo{t)) (38) 

and Xao{t) solves the SH system. All solutions of the SH system tend to the 
same limiting value X^. Thus by passing to a subsequence once more it can be 
seen that each sequence {t„} has a subsequence along which X{t) converges to 
X^. Thus in fact X{t) X,. 

If a solution of (|17p - (fT8)) does not converge to a stationary solution then it 
must converge to a periodic solution. Moreover, there are solutions where this 
happens whenever the parameters are such that the stationary point (x,, y,) is 
unstable. In this case it is more difficult to understand the behaviour of z and 
of the modified version of the SH model with variable c and C. 

5 Conclusions and outlook 

In this paper the properties of some mathematical models of parts of the NFAT 
signalling pathway have been investigated. One of the main results concerns a 
model describing the phosphorylation states of the transcription factor NFAT 
and the transport of these between the cytosol and the nucleus. Using chemical 
reaction network theory it was shown that for any value of the parameters in 
this system every solution converges to a stationary solution as t — > cxi and that 
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this stationary solution is uniquely determined by the total amount of NFAT in 
the cell. Moreover it was shown that in the case where transport between the 
two compartments is slow compared to the reactions within each compartment 
the stationary solution can be approximated by explicit expressions derived in 
|22j . It was exhibited in which way this model can behave as a switch. 

The parameters in the model describing the phosphorylation states include 
some which reflect the concentration of active calcineurin in the cytosol. It is a 
priori not clear that it is sufficient to assume that this concentration is constant 
and for this reason a model allowing a time-dependent calcineurin concentration 
was examined. The concentration of active calcineurin depends on the calcium 
concentration in the cytosol and this concentration is affected when a T cell 
is activated. The core model studied is a two-dimensional dynamical system 
closely related to a model for calcium dynamics introduced in [24]. Solutions 
of this model become periodic at late times for some values of the parameters 
and converge to a constant at late times for others. Criteria were obtained for 
the parameter values leading to these two different outcomes. It was shown 
that when the solution converges to a constant at late times this leads to a 
calcineurin concentration which also becomes constant and that the resultant 
densities of the phosphorylation states of NFAT in these cases are as in the case 
of constant calcineurin concentrations. 

In the course of this work the following additional questions arose. The 
criteria obtained for when exactly the system (IT71) - P^ does or does not have 
periodic solutions are far from exhaustive. In particular, one question left open 
is if there are parameter values for which the stationary solution is stable but 
there nevertheless exist periodic solutions. It was shown how the SH model is 
affected by an asymptotically constant input. It would be interesting to know 
how it is affected by an asymptotically periodic input. More generally it may 
be asked what can be said about a possible generalization of CRNT where the 
reaction constants are replaced by periodic functions. 

It was remarked that the results obtained for a model like the SH model, 
where mass action kinetics are used, may change essentially if another type 
of kinetics, such as Michaelis-Menten via mass action, is used. What happens 
when this is replaced by effective Michaelis-Menten kinetics? Here we are talking 
about different ways of modelling the same system of chemical reactions and it 
should be possible to relate them in a mathematically rigorous way. 

Ultimately, in understanding T cell activation, the interaction of the NFAT 
signalling pathway with the signalling pathways leading to other important tran- 
scription factors should be understood. Is it enough, when examining the cou- 
pled network, to confine attention to stationary solutions? Or so more com- 
plicated dynamical phenomena play a role? That they may do so is suggested 
by the results of ^ where it is shown that for instance NFAT and NFkB react 
to time-dependent calcium signals in different ways. A first step towards doing 
this could be to examine the influence of the dynamical behaviour within each 
of the individual pathways. One dynamical system describing NFkB signalling 
is given in [10]. The signal pathway leading from the T cell receptor to AP-1 
passes through a MAP-kinase cascade. Mathematical modelling of this pathway 
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revealing switch-like behaviour has been carried out in [T] . An insightful review 
of the possible contributions of theory and computation to the understanding 
of signalling pathways, with particular attention to T cells, has recently been 
given by Chakraborty and Das [3]. 

Acknowledgements I am grateful to Ria Baumgrass and Bernold Fiedler for 
helpful discussions. 

A Remarks on the model of Fisher et. al. 

The purpose of this appendix is to discuss the relation of the SH model discussed 
in the main part of the text to a model introduced by Fisher et. al. in jlOj. The 
building blocks in the latter model are NFAT and calcineurin. The NFAT may 
be unphosphorylated or fully phosphorylated and intermediate states are not 
included. The calcineurin may be active or inactive and the active form may 
bind to the different forms of NFAT. Concentrations of all these substances in 
the cytosol and the nucleus are considered and all of them may be transported 
between the two compartments. There are twelve concentrations in all. They 
satisfy a system of equations which, in neutral notation, can be written as: 

Vc 

kiX5 - k2Xi + kn — X2 - kig,xi + fcisxg - ki^xix^i 

Vn 

kiXQ — k2X2 + kis—xi — knX2 + ki^xio — kiQX2X4^ 

Vc 

Vc 

-kiix^x^ + ki2Xj + fcs — X4 - fce-Ts + ki^xg - kiex^xi + kiglxn - k2nxz 

Vn 



-kiiXeX4 + ki2Xs - ^5X4 + fcg 2:3 + ki^Xw - kiQXiX2 + kiglxi2 - k20X4 

Vc 

Vc 

-kiX5 + k2Xi - kiX^ + fca — xg - knx^x^ + ^12X7 

Vn 
Vn 

-kixe + k2X2 - fcsxe + k^ — X5 - knx^Xi + ^12X8 

Vc 

Vc 

kiiX^Xs. - ki2X7 + k'j X8 - fcsXy - /cisXy + kuXg 

Vn 



kiiXi^XA - ki2Xs + ^8 — X7 - krxs - ki^xs + kuxw 

Vc 



kl6XlX3 - ki5Xg + kijXj - ki^Xg + kg Xio ~ kiQXg 

Vn 



kiQX2XA - fcisXio + ^13X8 - ki^Xw + kig Xg - kgXiQ 

Vc 

Vc 

h — X12 - kexii - fcig/xii + fc2oa;3 

Vn 



-k5Xi2 + fee Xii - kiglxi2 + k20Xi 

Vc 



dxi 

IT 

dX2 

~dr 

dX3 

dt 

dxi 
IT 

dX5 

~dt 

dxe 

IT 

dxf 

~dt 

dxs 

~dt 

dxg 

~dt 
dxio 

dt 
dxii 

dt 
dxi2 

dt 
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Here all quantities other than the unknowns positive constants. To see 

the relation of these equations to those of [10] it suffices to note that they are 
written in the same order in both cases. There are conserved quantities 

Jl = Vn{xi +X5 + XT + Xq) + Vc{x2 + Xq + X^ + Xiq). (39) 

and 

■h = Vn{x3 + Xj + Xg + Xn) + Vc{x4 + Xs + XiQ + X12). (40) 

Jl is the total concentration of NFAT in all forms and J2 the total concentration 
of calcineurin. 

These equations do not appear to incorporate any switch-like mechanism 
and so it is difhcult to see how they could model the behaviour observed in [12] • 
The network structure contains rather little information about the nature of 
the substances involved and so it is to be expected that important information 
influencing the dynamics is encoded in the numerical values of the reaction 
constants ki. For instance, the fact that the central role of calcineurin in this 
context is that of a phosphatase is not visible from the equations. In [10] 
specific experimentally motivated values are chosen for the ki and the character 
of calcineurin as a phosphatase is encoded in the fact that fci3 is much greater 
than fci4. In view of this it is perhaps not surprising that CRNT is not too 
helpful in this case. There are sixteen complexes and only one linkage class. 
Since the rank of N cannot be more than the number of species, which is 
twelve, the deficiency must be at least three. Thus this system is far from the 
low deficiency situation where CRNT is most powerful. Nevertheless a couple 
of simple conclusions can be drawn. 

Because of the conserved quantities the reaction simplex is compact and all 
solutions exist globally in time. Suppose that a solution with Ji > and J2 > 
has an w-limit point on the boundary of the positive quadrant. Then arguing 
as in Section 2 we can conclude that there is a solution c* lying in the boundary 
and that at late times it has a constant number of non-zero components. Let 
rii be the number of non-zero components and uq — n ^ ni. Let Si and Sq be 
the corresponding subsets of S. It is impossible that there is a reaction whose 
left hand side is a species belonging to Si and whose right hand side includes a 
species belonging to Sq. Inspecting the pairs of species where there are reactions 
of this type in both directions shows that each of the following three sets are 
subsets of either So or Si 

S[ = {xi,x2,X5,xe}, S2 = {x3, 2:4, Xn, X12}, S'^ = {xr, xs, xg, xio}. (41) 

Moreover if S'^ C Si then Si — S. This possibility is ruled out by the fact that 
the solution lies on the boundary of the positive orthant. Inspecting the system 
once again and using the form of the nonlinear terms shows that if both S[ 
and S2 are subsets of Si a contradiction is obtained. Thus the only remaining 
possibilities are Si = 5( of Si — S2- However each of these implies the vanishing 
of one of the conserved quantities. It can be concluded that a solution for which 
Zi and Z2 are non-zero has no cj-limit points on the boundary. 



19 



It follows from the Brouwer fixed point theorerm as in Theorem 1.8.2 of [13] 
that there exists a stationary solution in the reaction simplex. Since it has been 
shown that this solution does not lie on the boundary it must lie in the interior. 
Thus this system has at least one positive stationary solution for any choice of 
the parameters and any biologically relevant initial data. On the other hand this 
argument gives no indication as to whether there might be multiple stationary 
solutions or more complicated dynamics such as periodic solutions. 
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